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The flow field created by swimming microorganisms not only enables their locomotion but also 
leads to advective transport of nutrients. In this paper we address analytically and computationally 
the link between unsteady feeding and unsteady swimming on a model microorganism, the spherical 
squirmer, actuating the fluid in a time-periodic manner. We start by performing asymptotic calcu- 
lations at low Peclet number (Pe) on the advection-diffusion problem for the nutrients. We show 
that the mean rate of feeding as well as its fluctuations in time depend only on the swimming modes 
of the squirmer up to order Pe 3//2 , even when no swimming occurs on average, while the influence of 
non-swimming modes comes in only at order Pe 2 . We also show that generically we expect a phase 
delay between feeding and swimming of l/8th of a period. Numerical computations for illustrative 
strokes at finite Pe confirm quantitatively our analytical results linking swimming and feeding. We 
finally derive, and use, an adjoint-based optimization algorithm to determine the optimal unsteady 
strokes maximizing feeding rate for a fixed energy budget. The overall optimal feeder is always the 
optimal steady swimmer. Within the set of time-periodic strokes, the optimal feeding strokes are 
found to be equivalent to those optimizing periodic swimming for all values of the Peclet number, 
and correspond to a regularization of the overall steady optimal. 



I. INTRODUCTION 



In order to be able to swim in viscous fluids, micro-organisms must undergo non-time-reversible sequences of shape 
changes referred to as swimming strokes (20l . [23l |32| . Through the no-slip boundary condition, these strokes induce a 
net flow field around the organism and a distribution of viscous stresses which lead to locomotion. This swimming- 
induced flow also impacts hydrodynamic interactions with neighboring organisms [ill . |29| or material boundaries 
[1, [l!| , the overall dynamics of suspensions of cells [H [H [M S| HI and the feeding ability of organisms [8L K35J. 

Cellular motility is essential to many biological functions, from reproduction [38| to escaping agressions 2, [l3|. It 
also allows organisms to travel toward better local environments for example to seek (or escape) light, nutrient, or 
heat. The performance of the particular stroke displayed by a single micro-organism, or that of a suspension of such 
swimmers, also results in the modification of the bulk stress and effective viscosity of a flow 0, HH, or of its mixing 
properties [jel [2ll . |34| , an effect that is suspected to play an important role on large-scale bio- mixing in the ocean for 
example |10j . 

The metabolism of many microorganisms relies on the absorption of diffusing nutrients present in their vicinity, 
ranging from dissolved gases and low-weight proteins, to more complex molecular compounds and, in the case of 
large organisms such as the protozoon Paramecium, smaller bacteria whose run-and-tumble motion is equivalent to 
a diffusive process at the scale of Paramecium (3|. For a particular microorganism, the impact of the stroke on its 
feeding ability can be thought of as twofold: (i) through the motility resulting from the stroke, the organism can travel 
toward nutrient-rich regions; (ii) by stirring nutrients in its immediate vicinity, the stroke-induced flow modifies, and 
possibly enhances, local concentration gradients. 

The competition of advective and diffusive effects on the dynamics of a particular nutrient is quantified in the Peclet 
number, Pe = Tdiff/V a d v , where raiff = a 2 / n and r ac i v = o,/U are the characteristic diffusive and advective time-scales 
respectively, where a, U and k are the typical size of the organism, the characteristic flow velocity, and the nutrient 
diffusivity, respectively. Depending on the nutrient considered, Pe can vary by several orders of magnitude, even for 
a given microorganism. 

Performing its stroke represents an energetic cost for the organism, as it must work against the fluid to overcome 
viscous dissipation. How far it can swim or how much nutrient it can absorb is therefore, in theory, limited by the 
finite amount of energy it has available. Considering that energy losses other than hydrodynamic can be accounted 
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for by a fixed metabolic efficiency, the optimization of the swimming stroke to maximize either motility or feeding can 
therefore be formulated as follows: for a fixed amount of energy available to deform its shape, what is the optimal 
stroke of a particular micro-organism maximizing either (i) the net displacement (optimal swimming problem) or (ii) 
the amount of a particular nutrient absorbed at the surface of the organism (optimal feeding problem)? In the latter 
case, the optimal stroke does not necessarily require a net displacement of the cell, as the organism can potentially 
just sit in a given location and stir the fluid around it. The optimal feeding stroke may also depend on the particular 
nutrient considered and the relative importance of advection and diffusion through the value of Pc. 

The optimization problems described above are closely linked to the question of optimality with respect to a specific 
biological function, which can take two different forms: optimal shape or optimal gait. In the former, one is interested 
in the optimal morphology of the swimmer (e.g. its aspect ratio, the use of flagella vs. cilia,...) and compares different 
species of microorganisms. In the latter, the focus is plac ed on a given organism, and the goal is to determine the 
sequence of body deformations that performs best [28ll30l l37l. I3ii - l4l| . 

In this work, we focus on the optimal gait of a particular swimmer model, the so-called squirmer. This canonical 
model, consisting of a spherical microorganism imposing a tangential velocity at its surface, was introduced as a 
so-called enveloppe model for ciliated microorganisms [If [HI- Ciliates, such as Paramecium, swim in viscous flows 
using the coordinated beating of a large number of small cilia distributed over their surface 0, 0] ■ I n the squirmer 
model, the flow field can be determined analytically through the projection of the stroke on orthogonal squirming 
modes. Because of its simplicity, this model has been used to study a large variety of problems related to swimming 
microorganisms, including hydrodynamic interactions [l5j . mixing (25| , suspension rheology [bH ]. collective dynamics 
and instabilities [HHl], and feeding [hHHES- 

Recently, Michelin & Lauga [28j determined the optimal time-periodic swimming strokes (i.e. those maximizing 
the swimming velocity for fixed energetic cost) of such a model microorganism, and identified their main properties. 
In a subsequent contribution, Michelin & Lauga (30| considered the optimization of the stroke for feeding in the 
particular case of a steady surface velocity. Although such strokes correspond to non-periodic displacements of the 
surface, the results shed some light on the link between swimming and feeding, and in particular it was shown that 
optimal swimming strokes and optimal feeding strokes were essentially identical regardless of Pe, a result that is not 
a priori intuitive due to the fundamental differences in the impact of swimming on feeding at low or high Pe: at low 
Pe, swimming only impacts marginally the nutrient distribution, but enables the organism to travel toward regions 
with richer nutrient content, while at high Pe, swimming also impacts feeding through stirring and strong advection 
of the nutrient in the vicinity of the organism surface. 

The validity of these conclusions, and in particular the intimate relationship between optimal swim ming and optimal 
feeding, remains however to be addressed in the general case of unsteady strokes. Magar & Pcdley |27| showed that 
in the particular limit of large Pe and small surface displacement, an equivalent steady problem could be defined. 
However, the unsteady effects of advection and diffusion in the general case of both finite swimmer displacement and 
finite Pe number remain unclear. In this paper, we specifically focus on the unsteady swimming problem. We first 
address analytically and computationally the link between unsteady feeding and unsteady swimming. We then derive, 
and use, an adjoint-based optimization algorithm to determine the optimal unsteady strokes maximizing feeding rate 
for a fixed energy budget. 

The paper is organized as follows. In £fl"fl the squirmer model is briefly presented, and the swimming and feeding 
problems are posed mathematically. In £]III[ the unsteady feeding rate is determined in the asymptotic limit of small 
Pe. The impact of the swimming stroke and of the Pcclct number on the feeding rate is further analyzed in £)IVI using 
numerical simulations, providing an important insight on the link between swimming and feeding. Section fVl presents 
the result of the stroke optimization with respect to feeding and conclusions and perspectives are finally presented in 

EH 

II. SWIMMING AND FEEDING OF A MODEL CILIATE 
A. The squirmer model 

The present work focuses on a particular model micro-organism, the squirmer, illustrated in figure [TJ It is a 
spherical organism of radius a which prescribes periodic tangential deformations of its surface S with a frequency 
ui, in order to swim in a viscous fluid of dynamic viscosity /i/ and density pf. The present analysis is restricted to 
purely axisymmetric deformations of S so that the swimming velocity is parallel to the axis of symmetry e x , with no 
rotation. In this paper, we will seek optimal strokes maximizing the feeding rate of the organism for a given amount 
of energy available during each period to perform its surface deformation (and possibly its swimming). This average 
rate of energy consumption, , is identified with the rate of work applied on the fluid by the swimmer at its surface, 
or, equivalcntly, the total mechanical energy dissipated in the fluid through viscous effects during one period. It is 
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FIG. 1: Swimming and feeding of a squirmer. A purely axisymmetric tangential velocity and a purely absorbing boundary 
condition for the nutrient are imposed at the surface of the swimmer. All variables are non-dimensional. 



related to the typical surface velocity scale % by 



?/ = J (I) 

The squirmer is swimming in a continuous suspension of a given nutrient (e.g. bacteria, large proteins/molecules, 
heat...) characterized by a far-field concentration Cx, and a diffusivity re, and advectcd by the flow created by the 
surface stroke. On the swimmer boundary, the nutrient is instantaneously absorbed and processed at the surface 
so that C = Cb, with Cb the equilibrium concentration at the surface determined by the processing mechanism. A 



more realistic, but more complex, boundary condition was proposed by Magar et at. |26[ and Magar & Pedley |27j . 
taking into account such effects as the resistance of the membrane to nutrient absorption, and the finite diffusion and 
processing time of the nutrient within the cell. The instantaneous nutrient uptake by the organism through diffusion 
at its boundary, <!>(£), is given by 



*(*) = y s re-dS. 



(2) 



In the case of a purely rigid sphere, with no advection, a steady nutrient flux is achieved through diffusion 4>o = 
47rare(Coo — Cb). In the following, we focus on the modification of the concentration field by the organism and define 
the rescaled concentration field c = (Coo — C')/(C 00 — C&). 

Three distinct time-scales are present in the problem: (i) a diffusive time-scale = a 2 /re, (ii) an advective time- 
scale r a = a/W and (iii) the stroke period t u = 2tt/uj, while only the latter two were present in the purely swimming 
problem |28l ] and only the first two in the steady feeding problem (30| . The Peclet number, Pc = Td/T a , is a measure 
of the relative importance of advective and diffusive effects near the surface of the squirmer, and is equal to 




Pc= — = ~\h- n (3) 



A second independent time-scale ratio can be defined either as a characteristic of the stroke, for example the relative 
velocity Ur = <fr/(aoj), or as a period- based Peclet number Pe w = a 2 w/re. In the following, all equations and 
quantities arc non-dimcnsionalizcd using a, uj, /if, and — Cb as reference quantities. 



1. Swimming problem 

Due to the small size of the organisms considered, the Reynolds number, Re = pftyfa/[if, a relative measure of 
inertia and viscous effects in the flow, is always much smaller than one, and the velocity and pressure fields satisfy 
Stokes equations. The swimming problem in the reference frame attached to the organism is therefore 



V 2 u = Vp, V ■ u = 0, 



(4) 



4 



with the boundary conditions on the swimmer surface and at infinity given by 

u = Ug (/i, t)eg at r = 1, (5) 
u — > — U(t)e x for r — > oo. (6) 

Note that the prescribed surface field, uf , is the stroke imposed by the organism and at the origin of both locomotion 
and stirring. The stroke is assumed to be axisymmctric, therefore the surface velocity only depends on [i = cos 6 and 
t, with 9 the polar angle measured from the swimming direction e x (figure [l}. In Stokes flow, the swimmer can not 
sustain any net hydrodynamic force, therefore we have 

/ [-pl+ (Vu + Vu T )] .nd5 = 0, (7) 
Js 

where n is the unit normal vector pointing into the fluid (n = e r here). Note that we have assumed the swimmer 
to be neutrally buoyant. The solution to the swimming problem in (|4j) — ([Tj) is obtained by decomposing the surface 
velocity onto the squirming modes 0, [28| 



(l*,t) = Y^a n (t)K n ( fJl ), (8) 



n=l 

with 

Kn(n) = -^—y/T^L'M, (9) 
n(n + 1) 

where L n {ij) is the n-th order Legendre polynomial. The values of the pressure field and strcamfunction are then 
obtained as 

oo 

p(r, fi,t) = poo + ^a n (t)P n (r, /i), (10) 

n=2 

oo 

ip(r, ^t) = Y^ a n (t)V n (r, n), (11) 

Tt = l 

with 

*n(r,M) = -^-(l-^)L'MMr), (13) 
n(n + 1) 

Mr) - Mr) = j(^^)forn> 2. (14) 

In the decomposition above, the first mode is the only one that contributes to the swimming motion (we have 
U(t) = a±(t) for all times) and is referred to as the swimming mode, or "treadmill". All remaining modes (including 
the so-called stresslet, n = 2, characterizing the modification of the bulk stress by the swimmer) correspond to higher 
order singularities in the far-field flow and do not contribute to the swimming motion. 
The dimensionless energetic cost, V, is computed as [28| 

P = i2^?=S^>' (15) 

with 

(2n + l) 2 

7l = 1 and 7 „ = ±- for n > 2, (16) 

Sn(n + 1) 

and is equal to the rate of working of the squirmer on the fluid through its boundary actuation or, equivalcntly, to 
the total energy loss through viscous dissipation in the fluid domain. In the following, we define (/) = ^ f(t)dt 
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as the time-averaging operator over one stroke period. With this definition, U = VV is the typical non-dimensional 
surface velocity of the swimmer. Following Lighthill [23T ]. the stroke swimming efficiency, r\ - or scaled energy cost -, 
is defined as the ratio of the energetic cost of pulling a rigid sphere with constant velocity (U) and the energetic cost 
of swimming at the same average velocity, obtained here as [281 ] : 

(U) 2 (a,) 2 

V = -Tfp- = —55 (17) 

n=l 



2. Feeding problem 

To evaluate the amount of nutrient absorbed at the surface of the organism, the non-dimensional advection-diffusion 
problem must be solved 

dc \ Pe 

— +u-Vc = V 2 c, wither —=, (18) 
at J VP 



together with the far-field behavior and purely absorbing boundary conditions on the swimmer surface (figure [1} 

c —> for r —> oo, (19) 
c = 1 for r = 1. (20) 

In equation l|18p the parameter e = ioa 2 jn can also be understood as the period-based Peclet number. The flow field, 
u, originates from the organism stroke and is obtained from the squirming mode amplitudes, ct n (t), using (fTTj) , (fT3")) 
and (fT4|) . The feeding performance of the stroke is evaluated using the ratio J(t) = $(t)/$o quantifying the net gain 
in nutrient uptake in comparison with the purely diffusive case (Pc = 0). The relative nutrient flux, J, is therefore 
non-dimensional and given by 



1 f 1 dc 



dfi. (21) 

r=l 



3. Eulerian vs. Lagrangian description 

A given periodic stroke, be it swimming or non-swimming, can be mathematically described following two different 
approaches: 

1. By prescribing at each instant, a periodic surface velocity on each point fixed in the swimmer frame, Ug(/j,,t), 
or equivalently a set of functions {a n (t)} n . We will refer to this description in the following as the Eulerian 
periodic stroke. 

2. By prescribing periodic trajectories, C(M0j0i °f material surface points labeled by their reference position on 
the sphere /io- We will refer to this description in the following as the Lagrangian periodic stroke. The surface 
velocity and mode amplitudes, a n (i), can then be obtained from £(no,t) as [28| 

1 d£ 

u e (a^o,t),t) = ^fat), (22) 

1 f 1 d 2 ^ 

a n(t) = ^ J L n [CCjUO.*)] g Mo g t d ^°- ( 23 ) 

In both descriptions, the flow velocity is periodic and completely determined by the periodic functions a n (t). 
However, in the Eulerian formulation, material surface points do not necessarily have periodic trajectories. Indeed, 
periodic Lagrangian strokes only represent a subset of periodic Eulerian strokes, namely the ones guaranteeing that 
every surface point comes back to its original position at the end of a full stroke period. Despite its shortcomings 
regarding the description of material point trajectories, the Eulerian approach has been the most popular for models 
of swimmers because of its simpli city , and in particular the possibility to consider steady strokes corresponding to 
steady surface and flow velocities [l(J [l2|, [ll| to cite only a few] . 
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B. Optimal swimming and optimal feeding 



For a given amount of energy available to perform a periodic stroke, an organism might have different optimal 
surface motions depending on the biological function of interest: migration (swimming problem) or nutrient uptake 
(feeding problem). A priori, those two objectives should lead to different optimal strokes, if anything because the 
optimal feeding stroke may depend on nutrient diffusivity through the value of Pe while the swimming problem does 
not depend on it. 

As emphasized earlier, a periodic stroke can be defined in two different ways, either from an Eulerian point of view 
(periodic flow field) or from a Lagrangian point of view (periodic material displacement). In our recent contributions, 
we presented the result of the optimal swimming problem (for both Eulerian and Lagrangian strokes) f28j and of 
the optimal feeding problem in the Eulerian steady framework only [30l ] . A brief summary of these results is first 
presented here. 

We start by remarking that, for the swimming problem, Eulerian optimal strokes are necessarily steady and each 
mode, a n , is independent of time. This is a direct consequence of the absence of history effect in the swimming 
problem: the swimming velocity and the energetic cost only depend on the instantaneous surface velocity. The 
optimal Eulerian stroke is then obtained by choosing the surface velocity distribution maximizing instantaneously the 
efficiency 77. From (|17[) we see that the Eulerian optimal swimming stroke is simply obtained by putting all the energy 
into the swimming mode, namely a n (t) = (L.i ■ The resulting treadmill swimmer, with an efficiency rj max = 50%, is 



In the case of the feeding problem, the presence of a time-derivative in the advection-diffusion equation introduces 
history effects, and the optimal Eulerian feeding stroke is therefore not necessarily steady. Focusing on the simplified 
problem of steady strokes, Michclin & Lauga [3 Of showed using adjoint-based optimization that the optimal steady 
feeding stroke is essentially the same as the optimal steady swimming stroke, a result which, surprisingly, remains 
true for all Peclet number. 

That result was not obvious a priori. The value of the mean feeding rate of the organism for a given stroke is a 
strong function of the diffusivity of the nutrient whose distribution around the organism is qualitatively different in 
the diffusive and advective regimes [2(|[3(|. The optimal feeding rate, (J) op t, depends strongly on Pe, but the stroke 
to achieve this optimal value does not. This result is important biologically as it implies that, for a given organism, a 
unique optimal stroke maximizes the nutrient uptake regardless of the details of its diffusive transport. For all Pe, and 
in the steady Eulerian framework, maximizing feeding and maximizing swimming are therefore equivalent problems. 

Although simpler conceptually and mathematically, the Eulerian framework is not appropriate to describe periodic 
deformations of a material surface, such as, for example, the strokes of ciliated cells. To impose periodicity of the 
surface motion, it is necessary to turn to the Lagrangian approach and to consider the unsteady swimming and 
feeding problems. Michclin & Lauga [28[ showed numerically that the optimal Lagrangian swimming stroke could 
be decomposed into two different parts: an effective stroke, dominated by the swimming mode, a\, and producing a 
forward velocity, and a recovery stroke during which material points (e.g. cilia tips) are brought back to their original 
position with front-like dynamics to minimize their (negative) impact on the swimming velocity. This front, or wave, 
is reminiscent of metachronal waves observed on ciliated organisms Q and results from a small phase-shift in the 
motion of neighboring surface points leading to global symmetry-breaking at the whole-organism level. When the 
squirmer model is used to represent a ciliate, the cilia length constrains the maximum displacement of the surface and 
therefore limits the ability of the swimmer to not only approach the optimal Eulerian stroke (treadmill) during the 
effective stroke but also to reduce the impact of the recovery stroke on the swimming motion. Using a constrained 
optimization algorithm, the direct relationship between swimming efficiency and surface displacement amplitude was 
obtained, and Michelin & Lauga [28[ showed that the optimal efficiency of 50% could be reached asymptotically. 

The optimization of the Lagrangian feeding stroke however remains at this point an open question; it is the focus 
of the present paper. The analysis of the nutrient uptake is first addressed analytically at small Pe. The general 
unsteady feeding problem is then considered numerically before turning to its optimization. 

III. UNSTEADY FEEDING AT LOW Pe: ASYMPTOTICS, SCALINGS, AND OPTIMUM 

In this section we focus on the feeding problem in the asymptotic limit of dominant diffusion (Pe <C 1). For a given 
stroke, this is equivalent to the asymptotic analysis of the advection-diffusion problem in the limit e = Pe/vP <C 1. 
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A. Steady and unsteady boundary layers 

For a steady velocity field, finding the asymptotic expansion of the scalar concentration, c, and surface flux, J, in 
the limit e< 1 corresponds to a variation on the classical mass transfer problem near a sedimenting sphere [T|, [2a . l30j . 
It is based on matching two different solutions for the scalar held c: near the surface of the sphere, diffusive effects arc 
dominant, and advection only appears as higher order corrections, while in the far-field, a balance of both advection 
and diffusion leads to the proper decay of c. 

In the case of an unsteady velocity field, both terms on the left hand-side of (fT8|) do not have the same scaling in 
the far-field. As a result the decay of the concentration field at infinity is not the same whether one considers the 
time-average of c or its fluctuations around the mean, and a double boundary layer problem must be considered: 

— in the near field, r = 0(1), diffusion dominates and the absorbing boundary condition (c = 1) at the surface of 
the swimmer is satisfied; 

— in the unsteady boundary layer (UBL), R = e 1 / 2 r = O(l), a balance between diffusive effects and rate of 
change of the local concentration ensures the proper far-field decay for the time-dependent fluctuations of the 
concentration field C(i?, p, t) = c(r, p, t); 

— in the steady boundary layer (SBL), p = er = 0(1), a balance between advection by the steady velocity field 
and diffusion ensures the far-field decay of the time-average concentration ^(p?^) = (c)( r :/')- 



B. Asymptotic problem formulation 

Decomposing the mode amplitudes, a n (t), as well as the concentration field, c, and feeding rate, J(t), into their 
Fourier components, we write 



a n(t) = ^2 CX niP Q lp \ C 



Cp(r,/i)e^, J(t)= -V 



ipt 



(24) 



The advection-diffusion equation becomes then 



oo oo 



in the near field, D • c p = e I ipc p 4- ot rLjq l n • c p - q J , 

\ n— 1 q— — oo / 



in the UBL, V ■ C p = ipC p + e 1/2 ^ ai, 9 A ■ C p - q + 0(e 



3/2 s 



in the SBL, @ ■ % = a lfi _£fi • % + 0(e 2 ). 

In (|25t - ([27l> . the following linear operators have been defined 



D 



1 



d ( n d 

7 a 
or 



d 



2n+ 1 



h = - 

In 



1 \ d 1 - A t 2 



/' 



dr r 
1 1 



d_ 

Of i 

i 



d 



2r 3 J dfj, ' 



di 



n-2\ d 



n(n+ 1) \r n + 3 r n+1 J dp: 



-p 



d_ _ (1-p 2 ) d 
OR R dp' 



(25) 

(26) 
(27) 

(28) 
(29) 
(30) 
(31) 



and T> (resp. 2!) is identical to D in (|28|) after replacing r by R (resp. p), and ££\ is defined as L\ after replacing R 
by p. The following boundary conditions must also be satisfied: 



\/p, Cp(r = 1) = Spfl, 
Vp ^ 0, C P {R oo) = 0, 
%(p-^oo) = 0. 



(32) 
(33) 
(34) 
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C. Matched Asymptotic Expansion 

A regular series expansion in e 1 / 2 of c p , C p and % is then performed up to 0(e 3 / 2 ). We write 

3 3 3 

Cp = E e9/2c P + °( £2 )' C p = ^ £ ^ + 0( £ 2 ), % = E e?/ X 9 + 0(e 2 ). (35) 

g=0 9=0 9=0 

At each order, c p and C p are to be matched for r — > oo and i? — > 0, while Co and & re to be matched in the limit 
R — > oo and p — > 0. 

Here, the non-homogeneous forcing (|32j) only acts on the steady-state component of the concentration field, and is 
transmitted to the time-dependent components by advection. Therefore, from the scalings of the different terms in 
equations ([25 ]) - ([26 ]) . 

Vp ^ 0, c p = O(ec ) and C p = O(e 1/2 C ). (36) 
1. Order O(l) 

At this order, advection is neglected and the solution is simply the steady diffusive solution C® = S Pi o/r, which 
satisfies both near-field and far-field boundary conditions. Therefore, C° = ^ = for all p. The resulting feeding 
rate is 

J P = S pfi + 0(e 1 / 2 ). (37) 
2. Order 0(e 1/2 ) 

Using (|36|) . = and C p = for all p ^ 0. The steady components cj, Cq, and % satisfy 

L> • Co = 0, (38) 
P • Cq 1 = 0, (39) 
=a 1 , Jf 1 -<. (40) 

Solving these equations and matching c p , C p and ^o up to C^e 1 / 2 ) leads to 

4 = 0, C p = ^, S*=0. (41) 
and the resulting feeding rate remains unmodified at this order. 

3. Order O(e) 

Next, the advection diffusion equation is expanded up to 0{e) in each region. 
— In the near field, r = O(l): 

D r 2 - V n 1 r°~ ^ p (\- M - V ( 2n + 1 K^ L "M f 1 _ 1 A f42 ^ 



n—1 x 7 n=2 

..2/ 



2r 2 

whose general solution satisfying the near- field boundary condition, c 2 (r = 1) = 0, is obtained as 



c 2 (r, ft) = a hp ft { A _ I _ + jr j n , p L n (fi) (-^ - rA 

\ ' n—1 ^ ' 



^ (2n+ l)a 7t|P L n (/x) / 1 | 1 2?i + 1 i 



n=1 



where j ntP are constants to be determined after matching with the UBL solution. 
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In the unsteady boundary layer, R = 0(1): 

-7-1 n 2 ■ n 2 r nl 

V-C p - \pC p = ai.pLi ■ C = —£2~, 
whose general solution compatible with the boundary condition at infinity (for p ^ 0) is 



r 2 _ ^o,o ,/ , ( 7i,o ,/ n ol\,o 



~ R 7o '° + M I i? 2 



C 2 = 2|P c -i ? v5 + Ai 



iai 



pi? 2 '' 7l 'f Vi? 2 ' R 
In the steady boundary layer, p = 0(1), the general solution of (|27[) is obtained as [l[ 



, for p 7^ 0. 



(44) 

(45) 
(46) 



»o = - ex P 
P 



g — m) ! 



(47) 



where the K 2 are constants to be determined in the matching process. Matching c p , C p , and ffg, up to O(e) leads to 



2 a i,o / 1 



(2n + l)an P L n ((i) 



s2 



E 

n=2 

Ql,p(l + /J.) 2 _ i Q i, P P 
2 ' P ~ pR 2 

1 / ai i0 (l + /u)p 
: P 6XP 1 2 



1 3 1 

2 4r 2 4r 3 

1 1 



2n + 1 



(rc +l)r™+ 2 nr n n(n + l)r n+1 J ' 



1- l + i?Vip e 



for p ^ 0, 



and the resulting feeding rate expansion is 



J P = S P ,o (l + ^)+0(, 3 / 2 ). 



(48) 

(49) 
(50) 

(51) 



Up to this order, we see that the results of the classical low-Pe asymptotic expansion for a steady velocity field are 
recovered and the mean feeding rate only depends on the average swimming velocity. In order to capture the leading 
order unsteady contribution to the feeding problem, the expansion must be carried out to the next order. 



4. Order 0(e 3/2 ) 

From ()21[) . we see that only the computation of the azimuthal average, c p , of the p-th Fourier component of the 



concentration field 

i r 1 

CpW = J 4(r,A*)dM, 

is necessary in order to compute the 0(e 3 ' 2 ) correction to the nutrient uptake. 

— In the near- field, taking the azimuthal average of (|25|) and using (|4ip . we have 

r 2 dr \ dr I 



(52) 



(53) 



whose general solution satisfying the boundary condition on the sphere is c p = a p (l — l/r), where a p is a constant 
to be determined by matching with the UBL solution. 
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In the unsteady boundary layer, taking the azimuthal average of (|26l) and using (|49[) , we get 

ip Cp = 2 ^ ai ' p ~~ q J l Cl ' 



q— — oo 



1 oo 

— £ «i,„ Ql ,,c-^. (54) 



3i? 

g— — oo 



This equation can be solved explicitly for Cp using the far-field boundary condition for the non-constant Fourier 
components and we get 



m— 1 

p R Z \pR 2y/i£ J ^ V 3mR I 



-RyJ i(p—m) — Ry/ir 



+ £ ^ — + " forp>l, (56) 

' \ 3ii / \ m p — m I 

m=p+l v 7 \ * / 

with defined for p < —1 using C_ p = C p . 

In the steady boundary layer, the equation for ^ 3 is identical to that at the previous order and the general 
solution takes the same form, see (l47l). 



By matching c p , C p , and ^0 U P to 0(e 3 / 2 ), the values of bo, a p , and a p can then be determined, and one obtains: 

(J) = 1 + + e 3/2 g ja^mp Q 2 
w 2 3 Jm 

m— 1 v 

J(t) - (J) = (-iVi) £ ^ + 0(e 2 ), , (58) 

p#0 



with 



p-i 

v m— 1 



E^ljm^l.m— p / 3/2 3/2 • / \3/2\ 



m>p+ 



(59) 



For a given stroke, the limit e <C 1 is equivalent to Pe <C 1 and the asymptotic expansion in terms of the Peclet 
number, Pe, can be obtained by substitution of e = Pe/vP in (|57| -([5^ |) . 

D. Discussion 

The asymptotic analysis obtained in (|57[) - (|59p provides some important physical insight into the relationship be- 
tween the swimming motion and the nutrient uptake on the surface of the swimmer. As for the steady case, the 
leading order advective correction to the feeding rate is linear in Pe and only depends on the average velocity of the 
organism P, HI]. At this order in Pe, there is a direct correlation between swimming and feeding and only the mean 
feeding rate is modified, fluctuations in time being negligible (higher order). 

The next order correction marks a fundamental difference between the steady and unsteady problems: in the steady 
case, all squirming modes contribute to the next correction at order Pe 2 [3(|. Instead, in the unsteady feeding problem, 
a new correction to J(t) (both its mean value in time and fluctuations) appears at order Pe 3 ^ 2 , which depends solely 
on the swimming velocity of the organism (through all the Fourier components, ai, m , of the swimming velocity, a±(t), 
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with no other squirming modes), and dominates the contribution of non-swimming modes that will only enter at 
order 0(Pe 2 ). For all time-periodic strokes, the instantaneous feeding rate is therefore completely determined up to 
(3(Pc 3 / 2 ) by the characteristics of the swimming velocity of the organism. 

This result has a major consequence for strokes that swim instantaneously (U(t) 7^ 0) but do not swim on average 
((£/) =0). In this case, the leading-order improvement to the feeding rate is solely governed by the zero- mean 
fluctuations of U(t). Non-swimming modes only contribute to higher order corrections, even if they have non-zero 
time averages. Consequently, for an organism that does not have a net swimming motion (e.g. a time-reversible 
swimmer), an instantaneous zero-mean swimming motion still presents a feeding advantage over stirring strokes 
where the cell stays in the same position at each instant (U(t) — 0). 

Our asymptotic expansion also provides some information on the relative phase of swimming and feeding. For an 
unsteady swimming velocity, U(t), with a single dominant Fourier component, the instantaneous feeding rate has a 
7r/4 delay on the swimming velocity (since — iVi = e -I7r / 4 in 155)) . A maximum in the feeding rate is therefore expected 
to take place after the peak swimming velocity, with a delay of l/8 th of a period. 

Note that the total nutrient flux is fully determined by the body velocity U(t) up to 0(Pe 3//2 ) . Whether the 
organism is swimming (force-free) or is an actuated rigid sphere (forced motion) does not actually come into play 
here. All the conclusions above are therefore valid for non-buoyant swimmers, but also for oscillating rigid spheres 
in Stokes flow, for which the present results represent a generalization of classical steady mass transfer results [l| to 
unsteady motions (see Appendix iBl for more details). 

In summary, our analytical results show that for low Pe, feeding is completely determined by swimming for any 
periodic stroke. Optimization of the feeding rate for a fixed amount of available energy is therefore equivalent in this 
limit to maximizing the swimming velocity under the same constraint, namely the swimming efficiency optimization 
problem. At low Peclet number, the Lagrangian optimal swimming and optimal feeding strokes are therefore identical, 
which confirms the result obtained in the steady framework by Michclin & Lauga [30| . In addition, similarly to the 
result for swimming, we get the result that at low Peclet number the optimal unsteady feeding problem is actually 
steady. This can be seen from (|57[) where the steady Fourier mode, 0:1,0, carries a higher weight than the other Fourier 
components compared to their relative importance in the rate of working. 



IV. UNSTEADY FEEDING AT FINITE Pe: SIMULATIONS 



To confirm the low-Pe results obtained analytically, we now turn to characterizing the feeding performance of 
different strokes for intermediate and large Pe. Eulcrian periodic strokes are determined by prescribing a n (t) for all 
n, while Lagrangian periodic strokes are described by giving the trajectories of material points 9 = &(6o,t) where 8 
is the current position of the material point and 6q its mean position. Alternatively, those strokes will be defined by 
fi = £(pL$,t), with /i = cos9. For illustration we consider three particular swimming and non-swimming Lagrangian 
periodic strokes: 



1. Stroke A is the numerical optimal swimmer identified in Michelin & Lauga 28| which has swimming efficiency 
f] sa 20%; 

2. Stroke B is a less efficient swimmer obtained using surface deformations in the form of a simple progressive 
wave: 

C(/i , t) = no + A(l - ill) cos(fc Mo - i), (60) 

with A = 1/3 and k = 1; 

3. Stroke C takes the same form as stroke B but with A = 1/3 and k = 0. Stroke C represents a time- reversible 
(or "reciprocal") deformation, and therefore has no net swimming motion, (U) =0. 

All three strokes display non-zero instantaneous swimming, but only strokes A and B show swimming on average. 
Stroke C differs thus from purely stirring strokes for which the organism is strictly still at each instant. The trajectories 
of material surface points are shown for strokes A, B and C in figured Mathematically, from the knowledge of £(/l«o, t), 
the mode amplitudes a n (t) are obtained using (|23[) . 
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FIG. 2: Trajectories, of surface material points for stroke A (left), stroke B (center), and stroke C (right). The correspond- 
ing time-averaged swimming velocity, (17), is equal to 0.33, 0.03, and respectively. The swimming and feeding performances 
of the three strokes are summarized in Table |T] 



A. Numerical solution of the advection diffusion problem 

For a given set of mode amplitudes, {a n (t)}, the advection-diffusion equation in (|18|) is solved spectrally in time 
for each azimuthal component of the concentration field 



oo oo 



p=0 



(61) 



k— — oo v— 



The functions c p (r) satisfy therefore the following systems of ordinary differential equations for p > and — oo < k < 



1 / d 2 2d p(p + 1) 



Pe Vdr 2 



with boundary conditions 



OO OO OO hi 

or ' 



X! X] "2 ( A mnp^n dr + B mnp ^ 



m—Q n—1 1— — OQ 



dr 



"p ( 
c k p {r 



c p {r = 1) = Sp t0 Sk,o, 
o) = 0. 



(62) 



(63) 
(64) 

are discretized on 
for more details], 



In (|6"2"j) . A mnp and B mnp are third-order scalar tensors defined in Appendix [5] Equations (|6"2")) -(1M 
an exponentially-stretched grid in r to concentrate points near the surface of the swimmer [see |3C 
and the solution { c p( r j)}^ kp - j i s then found iteratively. In typical simulations, the resolution used was AV = 120 
points for the r-grid, = 40-100 Legendre polynomials for the azimuthal dependence, N t = 16-128 points in time, 
and N a = 2-10 squirming modes to describe the swimming stroke. 

Alternatively, the advection-diffusion equation can be marched in time for each azimuthal component, c*(r, t), using 
an explicit time-stepping scheme for the advective terms and Crank Nicholson for the diffusion term. In the following, 
the advection-diffusion equation is solved spectrally in time except for strokes that do not swim on average (e.g. stroke 
C) for which the iterative algorithm does not converge properly or fast enough, and the time-marching approach is 
used in that case. 

Computationally, it is observed that the instantaneous nutrient flux converges rapidly with the number of squirming 
modes used to represent the swimming stroke, as shown in figure [3] The convergence is even faster for the average 
nutrient flux: describing stroke A with only the first two squirming modes significantly speeds up the computations 
while introducing an error smaller than 0.05% on the average feeding rate. Similar numerical tests performed on less 
efficient swimmers than stroke A (that is, swimming strokes for which mode 1 is not dominant) did not modify this 
observation significantly, and restricting the computation to only 2 or 3 squirming modes typically introduces an error 
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"'0 0.2 0.4 0.6 0.8 1 

t/T 

FIG. 3: Instantaneous nutrient flux at Pe = 5 for the optimal swimmer (stroke A) using an increasing number of squirming 
modes to numerically describe the stroke in the advection-diffusion solver: n a — 1 (solid), n a — 2 (dashed), n a — 4 (dotted) 
and n a = 8 (crosses). The error made on the average nutrient flux over a period is respectively 0.5%, 0.03%, 0.01% and 0.002%. 



Stroke 


(U) 


'/ 


(J) (Pe = 5) 


(J) (Pe = 30) 


A 


0.33 


22% 


1.97 


3.98 


B 


0.030 


1.3% 


1.33 


2.19 


C 





0% 


1.00 


0.99 



TABLE I: Swimming and feeding performance of strokes A, B, and C. 



smaller than 0.2%. This rapid convergence of the mean and fluctuating feeding rate is yet another indication that the 
swimming motion controls the feeding ability of the organism and higher-order modes only act as a small correction 
to the average feeding rate. 

B. Impact of the swimming stroke on the feeding performance 

Figures 2] and [5] show the concentration field around the squirmer for five successive and equispaced instants of a 
full period, for Pe = 5 (figure H]) and Pe = 30 (figure [5]), and for the three different strokes. For strokes A and B, at 
lower Peclet number, the nutrient concentration field only shows a weak front-back anisotropy as diffusion dominate 
over advection, confirming the observations on steady strokes of Magar et al. [2(| and Michelin & Lauga [3(| ■ As Pe 
is increased, sharper concentration gradients can be seen on the front of the squirmer. This results in an increased 
average feeding rate for increasing Pe as was observed for steady strokes [3(J. The main difference with the steady 
results is that in the unsteady scenario, the velocity of the squirmer changes (and possibly reverses sign) inducing a 
fluctuation in this front-back anisotropy and in the boundary layer thickness. For stroke C, which does not swim on 
average, the nutrient concentration field shows a strong isotropy, even at larger Pe, with much weaker concentration 
gradients resulting in a very weak modification of the nutrient uptake (J). 

Comparing the results obtained for the different strokes in Table HI we see that stroke A is clearly more efficient 
than strokes B and C from a feeding point of view, and stroke A also corresponds to a "better" swimmer. This is 
consistent with the increase of the feeding rate with the instantaneous swimming velocity that enables the formation 
of sharp concentration gradients in front of the squirmer. For stroke C, the periodic reversal of the swimming velocity 
over the period, and the absence of net displacement, results in the impossibility to maintain sharp concentration 
gradients at the front of the body and to swim toward regions with richer nutrient content, reducing its feeding ability 
significantly. 

Looking at the temporal variations of the swimming velocity and feeding rate throughout the stroke period (bottom 
frames of figures [4] and [5]) , a phase delay between the former and the latter is clearly identified for stroke A and B, 
and for all Pe considered. For stroke C, a similar delay is observed between the peaks in velocity magnitude (positive 
or negative) and the peaks in feeding rate: for this stroke, the feeding rate frequency is twice that of the swimming 




14 



Stroke A Stroke B Stroke C 




0.2 0.4 c 0.6 0.8 1 




0.5 1 0.5 1 0.5 1 



t/T t/T t/T 

FIG. 4: Top: Nutrient concentration around the organism at Pe = 5 for stroke A (left), stroke B (center) and stroke C (right). 
Bottom: Evolution in time of the feeding rate (solid) and swimming velocity (dashed). The dotted lines on the bottom figures 
indicate the time corresponding to each of the five top snapshots (ordered from top to bottom and left to right). 




FIG. 5: Same as figure g] with Pe = 30. 
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velocity because of the exact symmetry between the two half stroke periods. The presence of this time delay in all 
strokes is consistent with the results of the low-Pe asymptotic analysis in £11111 and can be interpreted as the time 
necessary for the concentration gradient (and possibly boundary layer) to reestablish at the front of the cell when its 
velocity starts increasing again. 

C. Impact of the Peclet number on the feeding performance 

It was observed previously that the value of Pe plays an important role in the feeding ability of the cell. This is 
investigated further here by looking at the impact of Pe on the instantaneous feeding rate for strokes A, B and C. The 
instantaneous feeding rate, J(t), is decomposed into its mean value, (J), the amplitude of its fluctuations in time, Ji, 
and its normalized profile, J(t), so we write 

J(t) = (J) + jJit), (65) 

where J\ = max(J) — min(J) and J(t) = (J(t) — (J))/J\. Similar quantities are also defined for the swimming 
velocity: (U), U\, and U. For a given stroke (A, B or C), the variation of these three quantities with Pe is displayed 
in figure [5] 

For swimming strokes, it is observed that, for low Pe, the modification in the mean feeding rate, (J) — 1, scales 
linearly with Pe (strokes A and B). This is consistent with the asymptotic analysis of Section ITU and with the steady 
results in Michelin & Lauga [30|. In such a diffusion-dominated regime, swimming enables the cell to sweep a region 
of fresher nutrients with an effective cross-section radius that is independent of the swimming velocity (because of the 
predominance of diffusion) and of the order of the size of the cell. At higher Pe, the reduced importance of diffusion 
over advection reduces the effective cross-section radius and (J) increases at a lower rate with Pe. For strokes with 
no net swimming motion (stroke C), the modification in the mean feeding rate scales as a higher power, Pe 3/2 , for 
Pe < 0.1, consistently with the results of the asymptotic analysis. 

For both swimming and non-swimming strokes, the amplitude of the feeding rate fluctuations, Jx, varies as Pe 3 / 2 
for Pe < 1, consistently with our asymptotic results. On figure [6] the fluctuations profile, J(t), is also represented and 
compared to the leading order prediction of the asymptotic analysis. We see a very good agreement at low Pe which 
persists even at high Pe for efficient swimming strokes such as stroke A. This confirms that the feeding rate (both 
its mean value and its fluctuations) is determined at leading order by the swimming mode and corrections from the 
other modes only play marginal roles. Again, a clear phase delay between the swimming velocity and feeding rate is 
observed for all Pe, and for the least efficient swimmers considered (B and C), this delay seems to increase with Pe. 

When Pe becomes large, another significant difference appears between strokes with zero (stroke C) or non-zero 
(strokes A and B) mean swimming velocity. For strokes A and B, the average feeding rate continues to increase with 
Pe, albeit more slowly. From the large-Pe steady results by Michelin & Lauga [3(|, we expect (J) to scale as Pe 1 / 2 , 
when the increase in feeding rate with swimming is driven by the concentration boundary layer thickness around the 
cell. In contrast, for non-swimming strokes, (J) reaches a maximum for a finite value of Pe (Pe c i ss 2) beyond which 
an increase in Pe actually results in a decrease of the feeding rate. Moreover, beyond a second critical value (Pe C 2 ~ 11 
for this particular stroke), the mean feeding rate falls below 1, and for large Pe, swimming actually penalizes feeding 
as it reduces the net feeding rate below the level of the purely diffusive regime (Pe = 0). This somehow surprising 
result can be understood as follows. In stroke C, the sphere swims forward during half of a period leaving behind it 
a nutrient-depleted wake. In the second half of the stroke, the cell swims backward into this region of poor nutrient 
concentration, resulting in a reduced flux at the boundary. 

D. The optimal unsteady stroke is steady 

As we discussed above, the optimal Eulerian swimming stroke is necessarily steady. The same conclusion can not 
be drawn a priori for the feeding problem due to the time-dependence of the advection diffusion equation (see £|II Bp . 
We saw however that it was true analytically at low Peclet number. Numerically, it also seems to hold as illustrated 
in figure [7J We performed numerical simulations on a large collection of unsteady Eulerian periodic and Lagrangian 
periodic strokes (8500 in total), ranging from very efficient to poor swimmers. For all values of Pe, the feeding rate 
is seen to be always less than that obtained with the optimal steady feeding stroke (treadmill). As for the optimal 
swimming stroke, the optimal Eulerian unsteady feeding stroke must therefore also be steady. Furthermore, figure [7| 
demonstrates that the more efficient the unsteady stroke is for swimming, the closer it can get to the optimal feeding 
rate. 




FIG. 6: Top: Dependence of the mean feeding rate, (J) (stars), with Pe and comparison with the asymptotic prediction in (|57[) 
(dashed). Center: Dependence of the peak-to-peak amplitude of the feeding rate fluctuations, J'(t) = J(t) — (J) (stars), with 
Pe and comparison with the asymptotic prediction in (|58p - (|59|) . Bottom: Rescaled (unit amplitude) feeding rate (solid) and 
velocity (dashed) time fluctuations; the asymptotic prediction for the feeding rate fluctuations at low Pe in (|58[) - l|59p is shown 
as a thick grey line. All results are plotted for stroke A (left), stroke B (centre), and stroke C (right). 



E. Feeding and swimming 



In the previous sections, a relationship between the swimming velocity and the feeding rate was clearly identified 
suggesting that at leading order, the mean feeding rate is determined by the swimming velocity and Pe. More precisely, 
and in the light of the steady results of Michelin & Lauga [30| , one expects the feeding rate to be determined by the 
swimming Peclet number, Pejj, defined as 

a(U) r- , . 

Pec; = = Pev/2^, (66) 



which measures the relative importance of advection of nutrients by the net displacement of the cell and diffusion. 
This is clearly the case at leading order for low Pe, as seen in ([57| . 
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Pe 

FIG. 7: Mean feeding rate as a function Pe for 8500 different swimming strokes (see text). The dashed line corresponds to the 
optimal steady feeding stroke (treadmill). For each stroke, the symbol color is related to its hydrodynamic efficiency, rj. 



In order to test the validity of this conjecture at higher Pe, we plot in figure [5] the mean feeding rate as a function 
of the "swimming Peclet number" Pec/, for the same large collection of unsteady strokes as in the previous section. 
All data points collapse rather well on a single curve, that corresponds exactly to the results for the steady treadmill 
swimmer [3(ij |. The agreement is particularly good for larger Pejj, corresponding to more efficient swimming strokes 
where the swimming motion dominates. The collapse of all the data points on that curve indicates that at leading 
order, for all strokes and all Pe, the mean feeding rate is determined by the mean swimming velocity. 

Figure [5] shows however that a significant number of points do not follow that leading order trend and are located 
above the grey treadmill curve. Indeed, for swimming strokes with poor efficiency (including those with Pejy = 0), the 
contribution from the mean swimming velocity to the mean feeding process is no longer dominant and the influence of 
other squirming modes, or from time- variations of the swimming velocity, cannot be neglected, so (J) remains strictly 
greater than one. 

V. OPTIMAL UNSTEADY FEEDING 

The results presented in the previous sections and in Michelin & Lauga f3fjj suggest that (i) swimming determines 
feeding, at least at leading order, and as a result (ii) optimal swimming and optimal feeding strokes are essentially 
identical. In this section, result (ii) is confirmed directly by performing an optimization of the swimming stroke 
maximizing the average nutrient uptake for a fixed energetic cost. The approach and methods presented below are 
based on the frameworks presented in Michelin & Lauga (H, H(| and generalized here to the unsteady feeding problem 
for periodic Lagrangian strokes. 

A. Adjoint optimization framework 

The rescaled nutrient concentration satisfies the advection-diffusion problem, (|18p - (j!20)) , and the mean feeding rate, 
(J), is given by 

( J ) = -(jZ I n-Vcds), (67) 
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where n = e r is the outward normal unit vector. Considering a small perturbation, <5u = ^2 n Sa n (t)u ( n >, in the 
velocity field, at leading order and for fixed Pe (or equivalently, fixed energetic cost) the resulting modification in 
mean feeding rate, (6 J) = S(J), is obtained at leading order as 

HJ) = -(IH(Sc)<1S), (68) 

where Sc is the resulting linear perturbation in the nutrient concentration field c satisfying 

d \ 5V 

— (Sc) + u • V<5c J - V 2 (5c = -eSu ■ Vc + — V 2 c, (69) 

with Dirichlct boundary conditions, 5c = 0, both on the surface of the swimmer and in the far-field. The last term 
in (|69[) guarantees that Pe = e\/V is constant and is obtained from Sa n and using (|15[) as 

5V = 2j2"/n{a n -5a n ). (70) 

n 

From (|69p , the change in mean feeding rate for constant Pe can be computed as 

S(J) =J2(an-Sa n ), (71) 

n 

where 

a„(t)= a ;(i)-^^((J)-K), (72) 
is the gradient of the feeding rate, at constant Pe, with respect to the n-th squirming mode amplitude and 

<^) = ^j a sVe-u^dO, H = (^J Vc-V.gdfi). (73) 
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In the previous equation, flf is the entire fluid domain, 

u («) i s 

the steady velocity field of the n-th squirming mode, 
and the adjoint field, g, satisfies the adjoint advection-diffusion problem 



together with boundary conditions 



§ + u.V. 9 



g — > 0, for r — > oo, 
g = 1, for r = 1. 



(74) 



(75) 
(76) 



A given Lagrangian periodic swimming stroke is defined by the trajectories of the surface material points, £(//o,i). 
The gradient of (J) with respect to the stroke, £(fj,o,t), is then the unique function F[^](fi , t) such that for any stroke 
perturbation the resulting modification in (J) is 



2tt Jo Jn, 



This gradient can be obtained directly from a n (t) as 

1 



F[£](m„t) = 



d 2 j 



d_ 

dt 



anLn djr 



(77) 



(78) 



and then projected onto the subspace of acceptable strokes (periodic trajectories, no displacement at the pole) [28j | . 
Note that although presented here for the particular case of a spherical swimmer, this optimization framework can 
easily be generalized to periodic swimming strokes of organisms with arbitrary shapes [30[ . 



B. Optimal feeding strokes 

Following Michelin & Lauga [28[ and in order to account for constraints on the stroke kinematics (introduced for 
example by a finite cilia- length-to-cell-size ratio), an additional constraint is included in the optimization algorithm 
to limit the maximum amplitude of angular displacements, ma x, of any surface point during the stroke. This 
optimization is performed using a steepest-ascent iterative optimization algorithm as described by Michelin & Lauga 
[30j . and the gradient of the feeding rate with respect to the swimming stroke is computed using the results from 
previous sections. 

Figure [5] shows the optimal strokes obtained for Pe = 5 and four increasing values of ma x- The optimal strokes 
consist in two different parts: an effective stroke where the surface of the squirmer stretches from front to back, 
enabling the swimming motion, followed by a recovery stroke where the material points (e.g. cilia tips) accumulated 
in the back side of the sphere are brought back to their original position with a front-like dynamics, reminiscent of the 
metachronal waves observed in ciliates. A wave velocity can be defined from the synchronization of the trajectories 
[28| . Notice also in figure [9] the phase delay between feeding and swimming predicted theoretically. Imposing tighter 
bounds on max results in a slower phase- velocity of the recovery stroke, in a smaller and steadier swimming velocity, 
and in a reduced efficiency (Table [TTJ) . This dichotomy of the optimal stroke and impact of the maximum displacement 
©max are essentially identical to that observed in the optimal Lagrangian swimming stroke by Michelin & Lauga p8j , 
for which it was observed that a continuous set of optimal strokes could be obtained for < © ma x < 90°, approaching 
asymptotically the optimal steady swimmer when max — > 90°. A similar behavior is observed on Figure fTUT a). 

The above conclusions are unchanged when performing the optimization at different values of the Peclet number, 
as shown in figure [TTJ For a given constraint on the maximum displacement max , the same strokes are obtained 
regardless of the value of Pe. These results confirm therefore that the optimal unsteady feeding stroke is essentially 
the same as the optimal swimming stroke, regardless of the value of the Peclet number. In both cases (swimming or 
feeding), the optimal Lagrangian stroke can be understood as a periodic approximation of the optimal steady stroke. 



VI. CONCLUSIONS 



In this paper we use asymptotic analysis and numerical computations to address the link between swimming and 
feeding for motile microorganisms. Using the mathematical model of spherical squirmers acting on the viscous fluid in 
a time-periodic manner, we first show analytically at low Pe that the mean rate of feeding as well as its fluctuations in 
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FIG. 9: Left: Lagrangian trajectories of four optimal feeding strokes, f(t), obtained for Pe = 5 and a maximum angular 
stretching of the surface equal to (a) max = 12°. (b) O m ax = 24°, (c) max = 35°, and (d) O max = 50°. Right: Time- variation 
for each of these optimal strokes of the instantaneous feeding rate, J(i) (solid), and swimming velocity, U(t) (dashed). The 
characteristics of these four strokes are summarized in Table [TTJ 



time depend only on the swimming modes of the squirmer up to order Pe 3 ^ 2 , even when no swimming occurs on average, 
while the influence of non-swimming modes come in later at order Pe 2 . We also demonstrate the existence of a phase 
delay between feeding and swimming of l/8th of a period. Using three illustrative stokes, we then employ numerical 
computations to confirm our asymptotic results and further demonstrate the relationship between swimming and 
feeding. Using adjoint-based optimization we finally determine numerically the optimal unsteady strokes maximizing 
feeding rate for a fixed energy budget. The overall optimal is always the steady swimmer. For time-periodic strokes, 
we find - as in the steady case - that the optimal feeding strokes are equivalent to those optimizing swimming - this 
result is true for all Peclet numbers even though the value of feeding rate strongly depends on the Pcclct number. 
As for the optimal unsteady swimming problem, optimal feeding strokes are therefore mathematical regularizations 
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(a) 
(b) 
(c) 
(d) 



"max V / 


(J) 


(U) 


n 


Pe 


e 


Pec/ 


12 


1.47 


0.048 


3.1% 


5 


26.2 


1.25 


24 


1.74 


0.141 


10% 


5 


15.9 


2.24 


35 


1.85 


0.221 


15% 


5 


12.3 


2.73 


50 


1.93 


0.295 


19% 


5 


10.5 


3.09 



TABLE II: Characteristics of the optimal feeding strokes obtained computationally for Pe 
displacements, © ma x, displayed on figure[9] 
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FIG. 10: Left: Maximum feeding rate, (J), as a function of the maximum angular displacement angle, max , in optimal strokes 
for Pe = 5. Right: Maximum feeding rate as a function of swimming efficiency, 77, for optimal feeding strokes obtained for 
Pe = 5 and various maximum displacement angles O max . The black stars in both figures correspond to the optimal steady 
stroke (treadmill). 



of the steady problem (treadmill) of overall maximum swimming and feeding performance. 

Clearly the problem studied here is idealized in many ways. The geometry is that of a sphere and the boundary 
conditions assume perfect nutrient absorption. These simplifications allow us however to develop a precise mathe- 
matical and computational description of the problem, both for the fluid and for the passive nutrient concentration. 
It is hoped that the biophysical insight developed in this study will be applicable to a wide range of problems in 
the realm of microorganism locomotion, e.g. in bacterial chemotaxis (at low Pe) or the feeding of plankton (at high 
Pe). One of the main modeling challenge for future work concerns the issue of shape changes. Most motile organisms 
display a Lagrangian deformation of their shapes. In this paper we have assumed that the deformations (the surface 
boundary conditions) always act tangentially to the organism surface, allowing the shape to remain that of a sphere. 
Clearly normal surface velocities would also need to be considered, and these are precisely the ones leading to changes 
in shape. The problem would then involve solving for the flow and nutrient concentration around a time-varying 
boundary. We hope that our study will inspire future work in this direction. 



Acknowledgments 

This work was supported in part by the US National Science Foundation through grant CBET-0746285. 



23 



(a)Pe= 0.5 



(6)Pe= 5 



(c) Pe= 25 
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FIG. 11: Left: Lagrangian trajectories of the optimal time-periodic feeding strokes obtained for max = 35° and (a) Pe = 0.5, 
(b) Pe = 5 and (c) Pe = 25. Right: Time-variation for each of these optimal strokes of the instantaneous feeding rate, J(t) 
(solid), and swimming velocity, U(t) (dashed). 



Appendix A: Definition of the -A mrL p and B mrL p tensors 

The coefficients A mnp and B mnp used in arc defined in terms of the Legendre polynomials as follow: 



A 
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2n(n + 1) 

They are easily computed using 

and the following recursive relations for n > 1 
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Appendix B: Unsteady heat /mass transfer around a sphere in Stokes flow 

In Section IIII1 the asymptotic expansion of the concentration distribution around a general squirmer and the 
resulting feeding rate, J(t), were obtained in the limit Pe < 1. The results obtained in (|58|) — (f5^|) also hold for any 
spherical object moving at velocity U(t), regardless of whether the sphere is swimming (zero net force) or a rigid 
sphere actuated by an external force, as we now show. 

Indeed, considering a generalization of the work of Acrivos <fc Taylor [l[ to unsteady particle velocity ai(t) = U(t), 
the velocity field around the sphere is given by the streamfunction 

, ax(t)(l - u 2 ) /3r , 1 A ,t,^ 
2 (j-^-Tr)- (B1) 

Following the same approach as in (|25j) takes the same form but (|29 |) — (f30 ]) become 

/31\<9l-/i 2 /31\ 
h = ~ [1 ~ Yr + 2^) Vfr ~ \ l Tr ~ 4^J ' (B2) 

l n = for all n > 2. (B3) 
In the same way, (|26p -(|2"T )l are slightly modified due to the contribution of the Stokeslet in the far-field: 

oo 
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OO 

+ e J2 a hA • + 0(e 3/2 ), (B4) 
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in the SBL, 9 ■ % = cti.o -Sfi ■ % + £"1,0 + 0(e 2 ), (B5) 

where C\ and ££\ remain unchanged from pip , and 

~ 3/i d 3(1 — fi 2 ) d 

Cl -2RdR + ^^^ (B6) 

and jSfi is obtained by replacing R by p in the previous equation. 

Following the approach of mill equations (j3"2"|) - ((5i?)) remain unchanged except: 



Equation (|42|) becomes 
Equation (|43|) becomes 
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These modifications do not impact the final result for the nutrient flux at the boundary. The expansion of the 
feeding rate for an oscillating sphere is therefore identical to that of the squirmer with same swimming velocity up 
to 0(Pc 3 / 2 ) in ([58]) - ((59)) . Looking at the corrections in the asymptotic expansion presented above, it appears that 
any far-field singularity in the velocity field (Stokcslct, etc..) will modify the near-field solution starting at O(Pe) 
and the unsteady boundary layer 0(Pe 3//2 ) but that such modifications will only affect the azimuthal fluctuations of 
the concentration and not its azimuthal average which determines the total feeding rate. Therefore, the asymptotic 
expansion of the feeding rate remains unchanged for any sphere moving at velocity c\i(t), regardless of the tangential 
velocity field applied on its surface, and regardless of the total force applied on the sphere. 

As a result, equations (|58[) (f59]) are a generalization to unsteady motions of the classical result on the heat and mass 
transfer on a sedimenting sphere [l[ , and the physical conclusions of £11111 are also valid in the case of a rigid sphere, 
in particular (i) the phase delay between the velocity and the mass transfer rate and (ii) an increase in mass/heat 
transfer scaling as Pe 3//2 for a sphere oscillating around a fixed mean position ((U) = 0). 
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